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Abstract 

The aim of the present paper is to efficiently describe the membrane 
potential dynamics of neural populations formed by species having a high 
density difference in specific brain areas. We propose a hybrid model 
whose main ingredients are a conductance-based model (ODE system) 
and its continuous counterpart (PDE system) obtained through a limit 
process in which the number of neurons confined in a bounded region of 
the brain is sent to infinity. Specifically, in the discrete model each cell 
of the low-density populations is individually described by a set of time- 
dependent variables, whereas in the continuum model the high-density 
populations are described as a whole by a small set of continuous vari¬ 
ables depending on space and time. Communications among populations, 
which translate into interactions among the discrete and the continuous 
models, are the essence of the hybrid model we present here. Such an 
approach has been validated reconstructing the ensemble activity of the 
granular layer network of the Cerebellum, leading to a computational cost 
reduction. The hybrid model reproduced interesting dynamics such as lo¬ 
cal microcircuit synchronization, travelling waves, center-surround and 
time-windowing. 
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1 Introduction 

Interesting phenomena in the brain often involve complex networks with an ex¬ 
tremely large number of neurons. The description at the microscopic level of 
the whole network, i.e., the modelling of each single neuron and synapse, would 
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lead to numerical models of prohibitive computational cost, even on the most 
advanced computers. The difficulties of such a description may be alleviated to 
some extent by identifying a hierarchy among interacting populations of neu¬ 
rons, and by using models with different resolution and cost for simulating the 
behaviour of different populations. Cell density may be a criterion to identify 
families of neurons and to partition the network in a multi-level manner, where 
each level corresponds to one or more species with comparable density. In the 
simplest situation of a two-level organization, this option leads to describe each 
neuron of the low-density population(s) by means of an ODE system, and to 
characterize the high-density population(s) by exploiting a PDE system that 
describes the family as a continuum. The hybrid model collects the ODE and 
the PDE systems, as well as the fundamental interactions among them. Several 
efforts, which have resulted in the formalisation of different models, have been 
made to understand and reproduce the activity of high-density populations by 
reducing the degrees of freedom from many, i.e., the variable states for each 
neuron, to few. Mean-field, neural mass and neural-field models are some of the 
results of various “passage to the continuum” approaches. A review concerning 
these models can be found in m-m The major difference between neural-field 
models - such as the one we are going to present - and the others lies in the 
fact that the formers account for the spatio-temporal evolution of the variables, 
rather than considering just the temporal evolution of them. A pillar formalisa¬ 
tion of a neural-field model is proposed in [3]-[35j-[36]. in which the macroscopic 
state variable is the mean firing rate. A more general neural-field model, not 
necessarily involving only firing rate variables, is presented in [32j . 

We obtain a continuum model for the action potential of a dense population 
of neurons by starting from a discrete model and letting the number of neurons 
tend to infinity while keeping them confined in a bounded region. We identify 
limit operators, acting on the continuous variables, describing specific interac¬ 
tions: in particular, electrical couplings (“gap junctions”) are modelled in the 
limit by the Laplace differential operator, as it has been rigorously justified in 
; on the contrary, chemical synaptic couplings produce non-local integral op¬ 
erators, i.e., spatial convolutions with suitable kernels (see e.g. Sect. 9.2 in m)- 
Once the expressions of both the discrete and the continuum model have been 
set, we describe in a fairly general form how the two models reciprocally interact, 
producing a hybrid model: aside of terms in the equations describing interac¬ 
tions between “homogeneous” (i.e., discrete-discrete, or continuous-continuous) 
variables, new terms are added to account for the “heterogeneous” interactions 
(i.e., between discrete and continuous, or continuous and discrete, variables). 

To validate our new method in a complete workflow we applied it to a re¬ 
alistic computational problem, the reconstruction of the Cerebellum granular 
layer network (GLN). This brain area shows a simple network structure yet 
capable of generating complex activity patterns. This network layer is densely 
populate by granule cells (GrCs) and sparsely by Golgi cells (GoCs) providing 
an optimal application for our modeling approach. The proposed hybrid model 
was specialized to the description of the interactions between GrC and GoC 
populations in the Cerebellum. Interesting dynamics such as local microcircuit 
synchronization, center-surround and time-windowing, as already described in 
a previous and more biologically detailed model [SO], are reproduced by the 
proposed model. Moreover, our model show the emergence of travelling waves 
of network activity elicited by a specific input configuration. 
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2 Materials and Methods 


2.1 The hybrid model 

In this section, in order to introduce the hybrid model, we first show how to 
model each individual neuron belonging to the same population. Here, intra¬ 
population communications are taken into account. Secondly, due to the fact 
that the number of neurons even in a small brain area is often huge, we perform 
a continuum limit of the discrete model that describes single neurons, obtaining 
a continuous model. Finally, we present the hybrid model in which the discrete 
and the continuous models interact with each other. 

Let us start by analysing how to describe the dynamics of each individual 
neuron i in the network, with i = 1, ■ ■ ■ N, where N is the number of neurons in 
a population. Precisely, we consider three variables: the voltage-like variable Vi, 
the recovery variable and the Si variable which describes the fraction of open 
channels in the synapses. In the most general fashion, each neuron is influenced 
by other neurons in the network by means of electrical and chemical synapses, 
and its dynamics is also driven by terms that describe the basic properties of 
neural excitability. All these ingredients are taken into account in the following 
general model: 

^=9{vi,ri), ( 2 . 1 ) 

ds^ 

at 

where, is the input current that accounts for electrical synapses, and 
is that for chemical synapses. In particular. 


■^gap = d ’ 

^syn 9syn,i ^ ^ WijSjiVi Usyn,j) j 


( 2 . 2 ) 


where Q(i) and B{i), resp., collect the indexes of neurons connected to the z-th 
one by means of electrical and chemical synapses, resp., Wij are positive weights 
describing the directed connection strength from j to z, d > 0 is the diffusion 
coefficient, gsyn,i > 0 is the synaptic efficacy, and Vgyn,] is the reversal poten¬ 
tial of the presynaptic neuron whose sign determines the synapse nature, either 
excitatory or inhibitory. In m and m, a detailed classification of synaptic 
reversal potentials, linked to distinct neurotransmitter/receptor pairs, is speci- 
hed. Furthermore, among the wide variety of models which describe the basic 
properties of neural excitability, we select the FitzHugh-Nagumo model [17] phe¬ 
nomenologically extracted from the biophysically-based Hodgkin-Huxley model. 
Thus, 

f{vi,ri) = -Vi{a - Vi){l - Vi) - n , 

9{vi,ri) = bvi - CTi . 


Here, a, b, c & R"*" are parameters chosen so that Vi is a fast variable and is a 
slow one. Finally, in the third equation in (2.1), a and /5 are positive parameters 
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describing the forward and backward rate constants for transmitters binding, vt 
is an o priori fixed threshold, and Hoo = Hoc (z) is the Heaviside function such 
that Hoo = 0 if z < 0 and Hoo = 1 otherwise. The model (2.1) obviously should 


be supplemented by suitable initial conditions for the variables {vi,ri,Si). 

In order to avoid prohibitive computational costs when the density of cells 
in a population is too high, we perform a “passage to the limit” as the number 
of neurons N tends to infinity in ( |2.1[ ). In this way, we capture the dynamics 
of a neuronal population as a whole by describing three continuous variables 


v{x,t), r{x,t) and s{x,t) (having the same meaning as in (2.1)), where x is the 
spatial variable. Specifically, in the limit case of —>• oo in a fixed and bounded 
spatial region H C M™, with m G {1, 2, 3}, the discrete model (2.1) leads to the 


following integro-differential system of equations (in which the t-dependence of 
each variable is ignored for simplifying notation): 


9syn w{x,y)s{y){v{x) - Vsyn{y))<iy 

Jn{x) 


^{x) = f{v{x),r{x)) + d*Av{x) 

^{x) = g{v{x),r{x)) 

f) c 

— (x) = a(l - s{x))Hooiv{x) - Vt) - /3s{x) , 
ot 

(2.4) 

supplemented by boundary conditions for v and initial conditions for v, r, s. 
Here, d* is the diffusion coefficient, gsyn > 0 is the synaptic efficacy, and TZ{x) 
denotes a region centered in x. The whole electrical synapse term, i.e. d*Av{x), 
is the result of two equivalent methods that lead to a non-trivial continuum 
limit, as shown in [^. On the other hand, the integral form of the chemical 

synapse term, i.e. g^yn f w{x, y)s{y){v{x) - 'Usyn(y))dy^ , is due to the fact 
that the set B{i) in (2.1) does not shrink to a point as N ^ oo, as explained in 


[T6] and [in- Furthermore, we refer to m for a discussion on the mathematical 
well-posedness of this model. Afterwards, in order to distinguish between the 
discrete and the continuum systems, variables in the continuum configuration 
will be indicated by Greek letters. 

As already mentioned in the introduction, by comparing the cell densities 
we may diversify the description of the populations in the network. Specifically, 
this comparison determines if a population is described by a set of discrete 
systems or by a continuous model. However, the key point is that neurons 
are linked to each others in a very intricate fashion depending on the brain 
areas. It follows that signal transmission among populations, in addition to 
intra-populations connectivity, is an important feature to be taken into account 
to explore the emergent network dynamics. The essence of the hybrid model 
lies in the interaction coupling terms among different populations. 

By considering for simplicity two populations only, on the one hand the set 
of cells in the low-density population is described by an ODE system: 


- 

= f{vi,ri) + ip{vi;Vj,Sj) + 4>(ni;w,cr) -(-, 
dr^ 

= g{vi,ri) , 

——- {Vi , 

dt 


(2.5) 
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where 


) d ^ ^ (^Vj Vi) 3syn ^ ^ ^synj) ^2 0^ 

jeQ{i) j&B{i) 

takes into account inputs from other cells belonging to the same low-density 
population, whereas 


^{vi;u,a) 


S^Loi^Xi^ Tsyn 



w{i, y)a{y){vi - u^yniy)) dy 


(2.7) 


describes the signal transmission coming from the continuous population. Here, 
Xi indicates the spatial position of the neuron labelled by i from the discrete 
family, whereas TZi is the region occupied by the neurons from the continuous 
family whose synapses influence neuron i. The term represents an external 
current coming from sources different from the two species here considered. On 
the other hand, the high-density population is characterized by a PDE system: 


9a; 

1h 

dp 

dt 

da 

~di 


F{uj,p) + V’(w,cr) -f T(a;;?;,s) -t-Xext , 
G{uj,p) , 

a(l - a)Hooiuj - lot) - Ida , 


where, similarly to (2.6), 


( 2 . 8 ) 


V'(i^,<^)(6 = - 7syn / w{C,y)a{y){oo{^) - u]syn{y))dy (2.9) 

Jn(0 

concerns interactions within the continuum population, while 


T(a;;u,s )(0 = d ^ {vj - uj{^)) - g^yn ^ w{(J)sj{uj{^) - Vsyn,j) 
jeQiO jeBiO 

( 2 . 10 ) 

describes the interactions between species, and Text = 2iext(0 is an external 
current. We call hybrid the model constituted by systems (|2.5|)-(|2.7|) and (|2.8|)- 

(2T0l). 


2.2 Application to the Cerebellum granular layer network 

The formalization of the hybrid model developed above is suitable for describing 
a variety of networks in the brain characterized by a large difference in their 
population densities. Among others, the olfactory bulb, the striatum, the gran¬ 
ular layers of the dorsal cochlear nucleus and the Cerebellum cortex are suitable 
to be efficiently represented with our new method. Out of these exempts the 
Cerebellum cortex is the most extensively studied and modelled network. 

The network structure can be abstracted following previous modelling works 
[ini Sni but keeping it sufficiently adeherent to the biological reality to show the 
versatility of the method to reproduce neural network dynamics observed in 
brain tissues. In the specific, the limit used to push to continuum the repre¬ 
sentation of the neuronal population with high density, could rise issues on the 


5 




reproducibility of network dynamics generated by a network composed by many 
small, independent units yielding unwanted diffusion of the activity across the 
network. A case that cannot yet be excluded in the nature of the real brain tis¬ 
sues but that was excluded in the biologically realistic simulations [SO]. More¬ 
over, the practical test should highlight the cooperation of the discrete units 
with the continuum model. 

Interests in Cerebellum dates back to the morphological studies carried out 
by Ramon y Cayal and Camillo Golgi, the electroencephalography studies car¬ 
ried out by Adrain [1] and the motor impairment manifest in World War I and 
II patients with cerebellar lesions studied by Holmes [TO]. Only later on, the 
hne Cerebellum structure inspired theories linking the network structure to a 
function starting, with Braitemberg [7], Marr [25], Albus [2] and Ito [20], a re¬ 
search work yet to be accomplished. Its peculiar structure comprehends series 
of highly regular, repeating units, each of which contains the same basic micro- 
circuit. The similarity in repeating units, from architectural and physiological 
perspectives, implies that different regions perform similar computational op¬ 
erations on different inputs. These inputs originate from different parts of the 
brain, spinal cord, and sensory system projecting directly into the Cerebellum. 
In turn, the Cerebellum projects to all motor systems. Despite the regularity 
of the Cerebellum facilitates its description, it remains a network able to gener¬ 
ate complex dynamics whose potentialities and functionalities are not yet fully 
understood. 

Few cellular populations in the Cerebellum cortex compose this geometri¬ 
cally regular network and are localised in three well distinct layers called molec¬ 
ular, Purkinje, and granular. The latter is densely populated by GrCs (density 
4.000.000/mm^) and sparsely by GoCs. The key point to support the appli¬ 
cation of our new modelling method is that the number of GoCs significantly 
differs from that of GrCs: GoCs are very few compared to GrCs [nnsniis] in 
the reason of about 1 : 400. Thus, by virtue of this strong density difference, the 
exploitation of combined discrete and continuum models becomes interesting. 
In particular, the variables {vi,ri,Si) describe each GoG through (2.5), while 
{uj,p,a) portray the GrG species as a whole by means of (2.8). We focused 
our test study to reproduce the transformation imposed to the input signals 
by the Gerebellum granular layer network (GLN). The ultimate output of the 
GLN provide excitatory input through their axons in the molecular layer to 
the Purkinje cells which constitute the only output pathway of the cerebellum 
cortex. The GLN is composed of two main network pathways, a feedforward 
path and a loop or feedback path, where both Granular cells (GrGs) and Golgi 
cells (GoCs) receive external excitatory inputs by the Mossy fibers (MFs) origi¬ 
nating from the precerebellar nuclei neurons. MFs excite both cell populations 
duplicating their input into two pathways. Along the feedback path MFs excite 
GrCs. These excite GoCs through the ascending axon and the parallel fibers 
(PFs), and GoCs, in turn, inhibit GrCs. In a compact writing: MF-GrCs-PFs- 
GoCs-GrCs. The second or feedforward path is constituted by the excitatory 
input from MFs to GoCs which terminates inhibiting GrCs. This pathway is 
MF-GoCs-GrCs. 

Inspired by assumptions in |28j and for modelling purposes, we consider the 
two populations belonging to two-dimensional parallel layers, as described in 
Fig. |2.1| The bottom one is constituted by the GrC continuum and the upper 
one collects GoCs. A third layer, above them, collects PFs. In reality, GoC 
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Figure 2.1: Connection topology between GrCs and GoCs from a postsynaptic neuron 
perspective: GrGs linked to the i-th GoC (left) and GoCs which are connected to the 
GrG at the point ^ (right). 


somata and GrCs are located in the just mentioned granular layer while the 
site where GoC dendrites receive input from the GrC axons (PFs) is in the 
molecular layer. 

Let us now define our model topology and connectivity in the GLN taking 
into account the fine structure of the biological network. The model was build 
to reproduce a GLN fraction with size 1500 [xm along the sagittal axis, 500 |xm 
on the transverse axis and 100 pm thick. However, in our representation the 
thickness of this flat volume is disregarded. In our model all spatial units are 
normalized to the network edge length. GoCs and GrCs are assumed to belong 
to two rectangular domains of size [0,3] x [0,1], one on top of the other (Fig. 
2.1). The projection of MFs inside the GLN shows an abundant parasagittal 
branching. Each MF innervates multiple cerebellum lobules. Within the lob¬ 
ule, local branching gives origin to small clusters of about 8 MF terminals in 
a rectangular area of 200 pm along the transverse axis and 150 pm along the 
sagittal axis [SUED]! data from the rat cerebellum. About 50 GrCs project their 
dendrites (maximum length 40 pm, mean length 13.6 pm) on a MF terminal. 
Therefore, the activation of a single MF should give rise to small spots of acti¬ 
vated GrCs with response intensity degrading from centre to periphery. In our 
model, the GrC population is represented as a continues sheet split into vertices 
by tessellation allowing the calculation of numerical solutions. In this config¬ 
uration, we assume that MF terminals provide excitatory input to a subset of 
the vertices. A diffusive term in the PDE spreads the input to the neighbour¬ 
ing vertices, the intensity fading to none for a distance equal to 40 pm. GoCs 
receive excitatory input from MF terminals from a wider area as GoC dendrites 
are longer than GrC dendrites and span a larger GLN volume |15j . 

Each GoC arborized axon reaches the granular layer throughout a paral¬ 
lelepiped volume [5] elongated along the sagittal direction, whose projection on 
the two-dimensional granular layer is a rectangle 650 pm long and 180 pm wide. 
A GoC sparsely inhibits GrCs lying inside the rectangle. GrC axons, i.e., PEs, 
ascend to the molecular layer, bifurcate, and run parallel to each other in ei¬ 
ther direction along the transversal axis, our x-axis, for a few mm crossing the 
GoC apical dendrites. Each PF synapses onto many GoC dendrites along its 
path. The GoC apical dendrites branch out in all directions sampling PF input 
from a cylinder in the ML represents in the original model by a circle of radius 
50 pm |15j . Therefore in our model, a GoC provide inhibitory input to all the 
GrCs located within a rectangle elongated along the sagittal axis, with length 
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1.3 and width 1/2.8. Each GrC influences all GoCs in a rectangle elongated 
along the transverse axis, covering the entire GLN extension, and narrow along 
the sagittal axis, covering 1/10 on either side of the PF wide stripe of the GLN 
(see Fig. 2.1). Notably, GoGs receive chemical excitatory synapses by GrGs. 


Furthermore, GoGs are linked among each other by gap junctions connecting 
their apical dendrites |33j . This electrical coupling is represented in our model 
by a diffusion term between the vertixes of the discrete model, i.e. in a first 
approximation a GoG is coupled only with its nearest neighbours. 

As already mentioned above, the Golgi cell system can be described by the 


model (2.5); the general expression of the functions and $, given in (2.6) and 


(2.7), takes here the following specific form: 


ip{vi;Vj,Sj) = d ^ (vj-Vi) 


E 


( 2 . 11 ) 


$(ui;a;,cr) =-7syn / w{i,y)a{y){vi 
■Ini 






Moreover, = /^ossy is excitatory input due to the MFs. Let us recall 
that, in (2.11), the reversal potential Wsyn may depend upon the presynaptic 


neurons and, thus, it must be included in the integral term. However, since only 
GrGs influence GoCs by means of excitatory chemical synapses, we suppose Wgyn 
to be constant and we bring it out of the integral, obtaining 

$(ui; w, cr) = -7syn ^ j w{i, y)a{y) dy^ (vi - uisyn) ■ 

The set TZi determines the area containing those GrGs which synapse onto the 
f-th Golgi cell. Taking into account that GrGs excite GoGs through the PFs, as 
specified above, we consider TZi as a thin rectangle whose horizontal symmetry 
axis is determined by the i-th cell projection (see Fig. |2.1[ left). The rectangle 
area is chosen by fixing a reasonably small parasagittal extension. 

Furthermore, concerning the coupling term between the two populations, it 
is well known that GrGs receive inhibitory chemical synapses from GoCs. Thus, 


the GrC continuum is described by the model (2.8), where the functions "0 and 
T, introduced in (|2.9[) and (|2.10|), take the following specific form: 


= 5Aa;(0 , 

T(a;;u,s)(0 = -5syn[ ^ 


(^(0 


( 2 . 12 ) 


As above, the reversal potential Usyn of presynaptic GoCs is supposed to be 
constant and then it is not involved in the summation. In order to consider 
inputs from Mossy Fibers, we set Text = 2imossy The discrete set collects 
the indexes of GoCs which influence the GrC continuum at the point thus 
describing the connection topology. According to [5] , a GoC axon reaches a rect¬ 
angular region in the granular layer, centered on its soma; therefore, a possible 
choice is: 

B{0 := {/ G N : x, G R^} , (2.13) 

where denotes such a rectangle centered on the projection of ^ on the GoC 
plane and oriented perpendicularly to the TZi direction (see Fig. 2.1, right). 











Since cells are described by the FitzHugh-Nagumo model, it is important to 
recall that the threshold is not involved in the single-neuron dynamics but it 
concerns presynaptic neurons at the synapse level. Indeed, when the presynaptic 
neuron exceeds the threshold, neurotransmitter release starts and influences the 
postsynpatic one. 

We close this section by a few words about the numerical treatment of our 
model. Concerning GoCs - which form a discrete set - they are placed at the ver¬ 
tices of a quasi-uniform triangulation of the upper domain; we use the triangular 
mesh generator BBTR, described in [4], with the mesh rehnement parameter 
chosen to yield 250 vertices (RehningOptions parameter set to 0.0035; Fig. 3.1). 
On the other hand, GrCs - which form a continuum in our model - are described 
by a set of partial differential equations that need to be discretized in space. To 
this end, we resort to a classical second-order centered finite difference method 
(see, e.g. m)- In particular, we consider 24000 nodes in the domain, lying 
on a regular grid, to represent the 300.000 GrGs. Therefore, using this grid 
size each vertex represents 12 or 13 GrGs. However, the results of the simu¬ 
lations turn out to be nearly independent on the GrGs grid rehnement, as it 
will be documented at the end of Section At last, time integration of the 
resulting coupled system of ordinary differential equations is accomplished by 
the MATLAB routine ODE45. We remark that the spatial discretization might 
be accomplished by hnite elements instead of hnite differences, thus allowing 
for the easy use of unstructured grids that adapt themselves to the formation 
of localized patterns; this will be object of future work. 


3 Results and Discussion 


3.1 Oscillatory activity in the granular layer 

Numerical simulations have been performed with the aim to validate the capa¬ 


bility of the hybrid network model, composed of (2.5)-t-(2.111 and (2.8)-t-(2.12), 
to reproduce the GLN activity simulated in a biologically realistic model [301. 
The network size is equivalent to a box with 500 p-m edge length along both the 
transverse and sagittal axes, and 100 pm thickness containing the cubic volume 
(500 pm edge length) of brain tissue simulated in [30] . 

Inspired by the orders of magnitude of the parameters in (^[30], we set: 


gsyn = 1, d = 0.05, 
for the Golgi cell discrete model, and 


P = 1^^°^ =01 

mossy mossy ^ 


(3.1) 


7sy„ = 0.05, ,5 = 0.005, X„,ossy(0=2:SoL = 0.1, 


(3.2) 


for the Granular cell continuous one. In particular, is applied to 3% of 

GrG nodes randomly chosen with uniform distribution. It is well known that 
MFs input GoGs, as well as GrGs. Since in the real GLN also GoGs receive 
excitatory input from MFs. We assume that 3% of GoGs receive /^oSy = 0.1. 
The current is applied for all t > 50 ms. In the meanwhile, MF current is 
maintained active to 3% of GrGs from f > 0 ms. Both GrGs and GoGs which 
receive the external current are randomly chosen with uniform distribution. 
The thresholds vt and ojt for GoGs and GrGs are both set to 0.5. The GoG 


9 








Figure 3.1; Domain decompositions obtained by exploiting the triangular mesh gen¬ 
erator BBTR in [4]. The RehningOptions parameter is set to 0.01, leading to a sparse 
mesh. 


potentials are described with vertical bars while GrC dynamics is shown with a 
continuous surface. 


A portrait of the GoC-GrC dynamics has been obtained by exploiting (2.5 H-( 2.11) 
and (2.8)+(2.12|, assuming the topology described by (2.13). The excitatory in¬ 
put delivered by MFs to GrCs drives their activity above threshold and induces 
an increase in GoC potentials. The subsequent inhibition elicited in GrCs by 
the GoC inhibitory feedback loop (MF-GrCs-PFs-GoCs-GrCs) suppresses the 
GrC activity and the cycle restarts. The same local microcircuit synchronous 
phenomena do arise in the biologically realistic model of reference |30j and it 
is a characteristic dynamics observed in the GLN in vivo [51] and models |22] . 

This dynamics is replicated with a specific period. Some significant snapshots 
are shown in Fig. |3.2| Moreover, at a later time of the simulation, t > 400 ms, 
the synchronous dynamics spontaneously converts in an interesting dynamics 
where excitatory waves travel in the whole domain involving both GoCs and 
GrCs, see Fig. |3.2| 


3.2 Center-surround and time-windowing 

Over the recent years several studies on the GoCs-GrCs network have been 
focused on the analysis of the integration of excitatory and inhibitory input by 
GrCs HU (TH] HU Uni ■ To further validate our modelling reconstruction we 
focused on reproducing the spatial and temporal interaction of excitation and 
inhibition in the GLN following the work presented in [50] • According to [IH1I3U] . 
the input delivered by a small bundle of MFs in the GLN elicits the activation 
of a cluster of GrCs, a spot 33 ± Sqm wide at 70% of the peak amplitude [25] . 
The spot is limited in size and in time by the properties of the feed-forward 
and feed-back inhibitory loops, due to the GoC integration properties and the 
arrangement of their axons. This phenomenon, defined center-surround and 
time-windowing in [12] . is the result of the mismatch between the small area 
excited by the MFs and the wider area inhibited by GoCs activated directly and 
indirectly, through GrCs, by the same MFs, in combination with the inherent 
delay of the inhibitory loops. 

This section is devoted to present noticeable center-surround and time- 
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Figure 3.2: Ensemble dynamics in the hybrid model. After an initial period of 
initialization (t < 180 ms), a synchronous phenomenon within each population arises 
and the network activity shows oscillations with a frequency of 13 Hz. After a few 
cycles (t > 350 ms) a travelling wave phenomenono arises. The oscillatory frequency is 
unaffected by the spontaneous emergence of the waves. The waves of network activity 
diffuse along the sagittal axis. GrCs are represented with the coloured continuous 
graph; GoCs are described with bars showing potentials multiplied by a factor 3 for 
graphical reasons. 
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Figure 3.3; Snapshots describing center surround phenomenon. No excitatory inputs 
from MFs reach GoCs. On the contrary, GoCs are excited by GrGs through the PFs. 
In turn, each active GoG inhibits GrGs lying on a thin rectangle. The stimulus was 
set on at t = 5 and set off at t = 15. 


windowing phenomena reproduced by models (2.51 and (2.8) 
significant comparisons with results in 


Furthermore, 


are presented. The aim of such com¬ 
parisons is to show that our hybrid model reproduces the same dynamics shown 
in reference articles in the field. Now we set: 


T 


= /GoC ^ Q 

mossy 


mossy 


(3.3) 


(3.4) 


5syn = 1, d = 0.005, 
for the Golgi cell discrete model, and 

7syn = 0.05, <5 = 0.5, X^oBsy(0=^So?sy = 0.2 

for the Granular cell continuous one. 

The activation of a spot in the network centre was achieved in the original 
model by activating the MF terminals located within a sphere of radius equal 
to 20 pm in the network centre. Gonsidering that the average length of GrC 
dendrites was set to 14 pm yielding an overall excited area of about 34 pm. 
In the simulations we run to reproduce the impulse response of the GLN, we 
mimicked this activation by providing excitatory input to GrG vertices within a 
circle with radius equal to 34 pm, 1/14.7 in our normalised units, and located in 
the network centre. The simulation reproduced an activated spot in the network 
centre of the same size shown in [30j (data not shown). In a second simulation, 
we increased the activated area to 50 pm in order to achieve a spot 33 ± 5 pm 
wide at 70% of the peak amplitude [23] as shown in Figure 3.3 at t = 12 ms. 


Let us assume that the connection topology is again described by (2.13) and 


let us consider MFs exciting GrGs in a circle, having radius 50, located in the 
centre of the domain. Figure 3.4 shows the GLN response to a stimulus set 
on at t = 5 and set off at f = 15 when the inhibitory connections were left 
active or blocked and their difference. The center-surround organisation of the 
inhibitory projections shapes the GLN response in space as it is evident from 
the enlargement of the active spot when those connections were blocked. 
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Figure 3.4: The GLN was activated by a 10 ms pulse delivered by MFs to GrCs. The 
GLN activation at time 7 ms from the initiation of the stimulus shows the maximal 
activation yielded by the excitatory input (E peak; upper left panel). After 6 ms 
the GLN activation fades off due to the emergence of the inhibitory feedback {E 2 
peak; upper right panel). After the block of inhibitory synapses the E 2 peak increses 
in amplitude and extension (7^2j]^; upper right panel). The amount of inhibition is 
calculated as the change in GLN amplitude due to the block of inhibition at 13 ms 
from the stimulus initiation (/; lower left panel). The center-surround is represented 
as in m as the difference from the E peak and the inhibition I (lower right panel). 
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In order to compare our result with those shown in Fig. 5 of |30j . let us 
reproduce in Fig. |3.4| the same computational steps to evaluate the effect of 
inhibition of the GLN activation. After the onset of the MF input the GLN 
initiates its response with 1 ms of delay reaching its maximal activation after 
2 ms, indicated as E peak in Fig. 3.4 After 6 ms the GLN activation fades 


off due to the emergence of the inhibitory feedback and we choose this time to 
measure the E 2 peak. Blocking the inhibitory synapses the E 2 peak increases 
in amplitude and extension (inhibition blocked: £’ 2 ^]^). As in |30j . the amount 
of inhibition I is calculated as the change in GLN activity amplitude due to 
the block of inhibition. The center-surround is represented as in |30j as the 
difference from the E peak and the inhibition /. 

Let us recall that our model constituted by (2.5) and (2.8) has been designed 
under strong simplifying assumptions that do not allow us to take into account 
the wide variety of phenomena in the single cell and in the whole network. 
Furthermore, the GrG layer has been described as a continuum. Nonetheless, the 
remarkable result obtained is that our model is able to reproduce the benchmark 
dynamics on the right-hand side, at least in the significant time range in which 
the center-surround phenomenon arises. Concurrently, the delayed activation 
of GoCs allows the response of GrCs to the stimuli to survive till the GoCs 
inhibition arises. This configures a time window where GrCs are allowed to 
transfer their activity to the subsequent network layers. The intervention of 
GoCs inhibition closes this window resetting the GrCs activity and making 
them ready to reliably transmit a new stimulus. 

Finally, we conclude the present section by stressing that the simulations 
provided in this paper turn out to be nearly independent on the GrCs con¬ 
tinuous population grid refinement. Indeed, focusing on the framework that 
describes the center-surround phenomenon, we exhibit a comparison among the 
solutions produced by the model with increasing number of nodes in the space 
discretization of the GrC population. In Fig. |3.5| the evolution in time of the 
membrane potential of two cells in the domain is shown, for different values of 
the spatial resolution. Convergence is clearly documented, thereby providing a 
sound background to the use of our numerical simulator. 


3.3 Computational comparison 

The computational performance of our new modeling method was assessed by 
running a simulation with an equivalent representation of a portion of the GLN 
in both simulators: NEURON nni and our hybrid model simulator. The simu¬ 
lation used as reference is the one reproducing the center-surround effect in [30] . 
On the one hand, in m the full model simulation of 250 ms of network activity 
(simulation run using the code available at |2^ required about 428 sec on a 
AppleMacBook Pro (Intel Core 2 Duo 2.93 GHz) for a network of 4001 GrCs 
and 27 GoCs. On the other hand, our network consists in 2116 GrC vertices 
and 27 GoCs. Considering an equivalent of 250 ms of activity, our simulation 
required about 71 sec. Therefore, our network simulator is roughly 6.7 times 
faster than the NEURON simulator. We must also recall that the output of our 
simulator is immediately available for visualization in MATLAB while the out¬ 
put generated by the NEURON simulator requires an additional 30 min of post 
processing to be visualized. This analysis quantitatively confirms the reduced 
computational cost of employing our simplified model instead of a detailed one, 
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Figure 3.5: Grid convergence for different structured grid resolutions of the continuous 
model for GrCs: 64, 256, 1024, 4096 nodes. The left panel shows the membrane 
potential lo plotted as a function of each time step for the node placed at (0.5, 0.5) in 
the centre of stimulation, i.e. subject to the excitatory input. The right panel shows 
the same plot for a node at (0.4, 0.4) outside the stimulated area and receiving only 
inhibitory input indirectly elicited by the feedback inhibitory loop. 


without losing information about such fundamental activity in time and space as 
the center-surronnd and the time-windowing. Let ns stress that improvements 
of our codes will lead to further time simulation savings. The most significant 
one will consist in translating our routines into a programming langnage that 
could be compiled rather than interpreted, i.e. C rather than Matlab, and 
in restructuring our code in order to take advantage of the multithreading or 
parallelization programming feature of the C programming language. 


4 Conclusions 


With the aim of efficiently describing the dynamics of nenronal populations hav¬ 
ing a strong density difference in specific brain areas, the present work collects 
new resnlts next to the ones presented in [9]. We started by stating the dis¬ 
crete conductance-based model (2.1) which describes the single cell membrane 
potential variation in time due to both electrical and chemical synapses. Af¬ 
terwards, the derivation of the continuum model was obtained. By letting the 


number of neurons tend to inhnity, we arrived at the complete model in (2.4). 


The two models, discrete and continuous, were then coupled to describe pop¬ 
ulations exhibiting in specific areas of the brain signfficant differences in their 
densities, allowing us to formalize the hybrid model. Specifically, each cell of 
the low-density population was modelled by the discrete model, whereas the 
whole high-density popnlation was described by the continnum model. Com- 
mnnications among populations, which translate into interactions among the 
discrete and the continuous models, are the essence of the hybrid model we pre¬ 
sented. Such an approach, which may lead to a significant computational cost 
redaction, was applied to the Golgi-Grannlar network in the Cerebellnm. Inter¬ 
esting dynamics such as synchronization, travelling waves, center-snrround and 
time-windowing were reproduced by the hybrid model. The two latter dynamics 
were compared with recent results in literature devoted to this specific network, 
confirming the capability of onr approach to reprodnce significant dynamics. 
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By proceeding on the path here traced, some improvements should be taken 
into account in a forthcoming work. A major objective should concern how much 
the network behaviours here reproduced are related to the specific properties 
of the FitzHugh-Nagumo single-cell description. Moreover, one should evaluate 
if a different single-cell model is able to reproduce other significant behaviours 
such as resonant dynamics. Finally, in order to make the model more adherent 
to the reality, a future work should include the plasticity in communication 
strength among neurons. 
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